1 Read data

headings_em <- read_csv("data/raw/headings-emilian.csv")
headings_eo <- read_csv("data/raw/headings-esperanto.csv")

emilian <- read_csv("data/raw/emilian-clean.csv", skip = 1, col_names = headings_em$new, na = c("", "NA", "na")) %>%
  mutate(
    understand = factor(understand, levels = c("NO", "AL", "50/50", "G", "VG")),
    speak = factor(speak, levels = c("NO", "AL", "50/50", "G", "VG"))
  )

esperanto <- read_csv("data/raw/esperanto-clean.csv", skip = 1, col_names = headings_eo$new, na = c("", "NA", "na")) %>%
  mutate(
    age = str_remove(age, "jaroj"),
    age = str_replace(age, "naskiĝis en la 1995a", "25"),
    age = as.numeric(age),
    understand = factor(understand, levels = c("NO", "AL", "50/50", "G", "VG")),
    speak = factor(speak, levels = c("NO", "AL", "50/50", "G", "VG"))
  )

2 Emilian

2.1 Participants

emilian %>%
  ggplot(aes(gender)) +
  geom_bar()

emilian %>%
  ggplot(aes(age)) +
  geom_histogram()

emilian %>%
  ggplot(aes(age)) +
  geom_density()

emilian %>%
  ggplot(aes(education)) +
  geom_bar()

emilian %>%
  ggplot(aes(age, education)) +
  geom_jitter(height = 0.2, alpha = 0.5)

emilian %>%
  count(profession) %>%
  ggplot(aes(reorder(profession, -n), n)) +
  geom_bar(stat = "identity")

emilian %>%
  count(languages_family) %>%
  ggplot(aes(reorder(languages_family, -n), n)) +
  geom_bar(stat = "identity")

emilian %>%
  count(languages_parents) %>%
  ggplot(aes(reorder(languages_parents, -n), n)) +
  geom_bar(stat = "identity")

2.2 Language

2.2.1 Understand

emilian %>%
  ggplot(aes(understand, fill = understand)) +
  geom_bar() +
  scale_fill_brewer(type = "div") +
  theme_dark()

emilian %>%
  ggplot(aes(understand, fill = gender)) +
  geom_bar()

emilian %>%
  ggplot(aes(understand, fill = gender)) +
  geom_bar(position = "fill")

emilian %>%
  ggplot(aes(age, fill = understand)) +
  geom_histogram(binwidth = 5) +
  facet_grid(understand ~ .)

emilian %>%
  ggplot(aes(understand, fill = profession)) +
  geom_bar()

emilian %>%
  ggplot(aes(understand, fill = profession)) +
  geom_bar(position = "fill")

2.2.2 Speak

emilian %>%
  ggplot(aes(speak, fill = speak)) +
  geom_bar() +
  scale_fill_brewer(type = "div")

emilian %>%
  ggplot(aes(speak, fill = gender)) +
  geom_bar()

emilian %>%
  ggplot(aes(speak, fill = gender)) +
  geom_bar(position = "fill")

esperanto %>%
  ggplot(aes(age, fill = speak)) +
  geom_histogram(binwidth = 5) +
  facet_grid(speak ~ .)

emilian %>%
  ggplot(aes(speak, fill = profession)) +
  geom_bar()

emilian %>%
  ggplot(aes(speak, fill = profession)) +
  geom_bar(position = "fill")

2.2.3 Read and write

emilian %>%
  ggplot(aes(read_write, fill = read_write)) +
  geom_bar()

emilian %>%
  ggplot(aes(read_write, fill = gender)) +
  geom_bar()

emilian %>%
  ggplot(aes(read_write, fill = gender)) +
  geom_bar(position = "fill")

emilian %>%
  drop_na(read_write) %>% 
  ggplot(aes(age, fill = read_write)) +
  geom_histogram(binwidth = 5) +
  facet_grid(read_write ~ .)

emilian %>%
  ggplot(aes(read_write, fill = profession)) +
  geom_bar()

emilian %>%
  ggplot(aes(read_write, fill = profession)) +
  geom_bar(position = "fill")

2.2.4 Attitude

emilian %>%
  select(educated:familiar) %>%
  pivot_longer(educated:familiar, names_to = "feature", values_to = "rating") %>%
  ggplot(aes(as.factor(rating), fill = as.factor(rating))) +
  geom_bar() +
  scale_fill_brewer() +
  facet_grid(. ~ feature)

emilian %>%
  select(educated:familiar) %>%
  pivot_longer(educated:familiar, names_to = "feature", values_to = "rating") %>%
  ggplot(aes(feature, fill = as.factor(rating))) +
  geom_bar(position = "fill") +
  scale_fill_brewer()

2.3 Urban vs Rural

emil_rur <- emilian %>%
  mutate(
    ru_ur = ifelse(
      str_detect(birth_place, "-RU"), "rural",
      ifelse(
        str_detect(birth_place, "-UR"), "urban",
        NA
      )
    )
  )

emil_rur_clean <- emil_rur %>%
  select(id, understand, speak, read_write, educated:familiar, ru_ur) %>%
  mutate(
    understand = ordered(understand, levels = c("NO", "AL", "50/50", "G", "VG")),
    speak = ordered(speak, levels = c("NO", "AL", "50/50", "G", "VG")),
    across(educated:familiar, ~ as.ordered(.x))
  ) %>%
 drop_na()
emil_rur %>%
  drop_na(ru_ur) %>%
  ggplot(aes(ru_ur, fill = understand)) +
  geom_bar()

emil_rur %>%
  drop_na(ru_ur) %>%
  ggplot(aes(ru_ur, fill = understand)) +
  geom_bar(position = "fill")

3 Esperanto

3.1 Participants

esperanto %>%
  ggplot(aes(gender)) +
  geom_bar()

esperanto %>%
  ggplot(aes(age)) +
  geom_histogram()

esperanto %>%
  ggplot(aes(age)) +
  geom_density()

esperanto %>%
  ggplot(aes(education)) +
  geom_bar()

esperanto %>%
  ggplot(aes(age, education)) +
  geom_jitter(height = 0.2, alpha = 0.5)

esperanto %>%
  count(profession) %>%
  ggplot(aes(reorder(profession, -n), n)) +
  geom_bar(stat = "identity")

esperanto %>%
  count(languages_family) %>%
  ggplot(aes(reorder(languages_family, -n), n)) +
  geom_bar(stat = "identity")

3.2 Language

3.2.1 Understand

esperanto %>%
  ggplot(aes(understand, fill = understand)) +
  geom_bar() +
  scale_fill_brewer(type = "div") +
  theme_dark()

esperanto %>%
  ggplot(aes(understand, fill = gender)) +
  geom_bar()

esperanto %>%
  ggplot(aes(understand, fill = gender)) +
  geom_bar(position = "fill")

esperanto %>%
  ggplot(aes(age, fill = understand)) +
  geom_histogram(binwidth = 5) +
  facet_grid(understand ~ .)

esperanto %>%
  ggplot(aes(understand, fill = profession)) +
  geom_bar()

esperanto %>%
  ggplot(aes(understand, fill = profession)) +
  geom_bar(position = "fill")

3.2.2 Speak

esperanto %>%
  ggplot(aes(speak, fill = speak)) +
  geom_bar() +
  scale_fill_brewer(type = "div")

esperanto %>%
  ggplot(aes(speak, fill = gender)) +
  geom_bar()

esperanto %>%
  ggplot(aes(speak, fill = gender)) +
  geom_bar(position = "fill")

esperanto %>%
  ggplot(aes(age, fill = speak)) +
  geom_histogram(binwidth = 5) +
  facet_grid(speak ~ .)

esperanto %>%
  ggplot(aes(speak, fill = profession)) +
  geom_bar()

esperanto %>%
  ggplot(aes(speak, fill = profession)) +
  geom_bar(position = "fill")

3.2.3 Read and write

esperanto %>%
  ggplot(aes(read_write, fill = read_write)) +
  geom_bar()

3.2.4 Attitude

esperanto %>%
  select(educated:familiar) %>%
  pivot_longer(educated:familiar, names_to = "feature", values_to = "rating") %>%
  drop_na() %>%
  ggplot(aes(as.factor(rating), fill = as.factor(rating))) +
  geom_bar() +
  scale_fill_brewer() +
  facet_grid(. ~ feature)

esperanto %>%
  select(educated:familiar) %>%
  pivot_longer(educated:familiar, names_to = "feature", values_to = "rating") %>%
  drop_na() %>%
  ggplot(aes(feature, fill = as.factor(rating))) +
  geom_bar(position = "fill") +
  scale_fill_brewer()

4 Locations

if (file.exists("./data/raw/geo.csv")) {
  cat("Reading geocoding...\n")
  geo <- read_csv("./data/raw/geo.csv")
} else {
  birth_em <- emilian %>% select(birth_place_it) %>% unique()
  geo <- geocode(birth_em, city = birth_place_it, method = "osm", verbose = TRUE)
  write_csv(geo, file = "./data/raw/geo.csv")
}
## Reading geocoding...
europe <- ne_countries(continent = "Europe", returnclass = "sf", scale = "medium")
ggplot() +
  geom_sf(data = europe) +
  geom_point(data = geo, aes(long, lat)) +
  coord_sf(xlim = c(7, 14), ylim = c(43, 47))

5 Modelling

lang_use <- bind_rows(
  emilian %>% mutate(language = "emilian"),
  esperanto %>% mutate(language = "esperanto")
) %>%
  select(id, age, understand, speak, read_write, language, educated:familiar) %>%
  mutate(
    understand = ordered(understand, levels = c("NO", "AL", "50/50", "G", "VG")),
    speak = ordered(speak, levels = c("NO", "AL", "50/50", "G", "VG")),
    across(educated:familiar, ~ as.ordered(.x)),
    mean_comp_raw = as.numeric(understand) +
      as.numeric(speak) +
      ifelse(read_write == "Yes", 5, 0),
    mean_comp = mean_comp_raw / 15
  ) %>%
 drop_na()

attitudes <- lang_use %>%
  select(educated:familiar)

5.1 MCA

attitudes_mca <- MCA(attitudes, graph = FALSE)
attitudes_dims <- attitudes_mca[["ind"]][["coord"]]
fviz_screeplot(attitudes_mca, addlabels = TRUE, ylim = c(0, 15))

fviz_mca_biplot(attitudes_mca, repel = TRUE)

fviz_mca_var(attitudes_mca, choice = "mca.cor", repel = TRUE)

fviz_mca_var(attitudes_mca, col.var = "black", shape.var = 15, repel = TRUE)

fviz_mca_var(attitudes_mca, col.var = "cos2", gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"), repel = TRUE)

fviz_cos2(attitudes_mca, choice = "var", axes = 1:2)

fviz_contrib(attitudes_mca, choice = "var", axes = 1, top = 15)

fviz_contrib(attitudes_mca, choice = "var", axes = 2, top = 15)

fviz_contrib(attitudes_mca, choice = "var", axes = 1:2, top = 15)

fviz_mca_ind(attitudes_mca, label = "none", habillage = "trustworthy")

fviz_mca_ind(attitudes_mca, label = "none", habillage = "friendly")

fviz_mca_ind(attitudes_mca, label = "none", habillage = "kind")

fviz_mca_ind(attitudes_mca, label = "none", habillage = "familiar")

5.2 Understand

lang_use$dim_1 <- -attitudes_dims[, 1]
get_prior(
  understand ~
    language * dim_1,
  data = lang_use,
  family = cumulative()
)
##                 prior     class                    coef group resp dpar nlpar
##                (flat)         b                                              
##                (flat)         b                   dim_1                      
##                (flat)         b       languageesperanto                      
##                (flat)         b languageesperanto:dim_1                      
##  student_t(3, 0, 2.5) Intercept                                              
##  student_t(3, 0, 2.5) Intercept                       1                      
##  student_t(3, 0, 2.5) Intercept                       2                      
##  student_t(3, 0, 2.5) Intercept                       3                      
##  student_t(3, 0, 2.5) Intercept                       4                      
##  bound       source
##             default
##        (vectorized)
##        (vectorized)
##        (vectorized)
##             default
##        (vectorized)
##        (vectorized)
##        (vectorized)
##        (vectorized)
understand_priors <- c(
  prior(normal(0, 3), class = Intercept),
  prior(normal(0, 1), class = b)
)

understand_pchek <- brm(
  understand ~
    language * dim_1,
  data = lang_use,
  family = cumulative(),
  prior = understand_priors,
  sample_prior = "only",
  backend = "cmdstanr",
  cores = 4,
  file = "./data/rds/understand_pchek"
)

conditional_effects(understand_pchek, effects = "dim_1", conditions = make_conditions(lang_use, "language"), categorical = TRUE)

understand_bm <- brm(
  understand ~
    language * dim_1,
  data = lang_use,
  family = cumulative(),
  prior = understand_priors,
  backend = "cmdstanr",
  cores = 4,
  file = "./data/rds/understand_bm"
)
conditional_effects(understand_bm, effects = "dim_1", categorical = TRUE, conditions = make_conditions(lang_use, "language"))

5.3 Speak

speak_bm <- brm(
  speak ~
    language * dim_1,
  data = lang_use,
  family = cumulative(),
  # Same priors as understand_bm
  prior = understand_priors,
  backend = "cmdstanr",
  cores = 4,
  file = "./data/rds/speak_bm"
)
conditional_effects(speak_bm, effects = "dim_1", categorical = TRUE, conditions = make_conditions(speak_bm, "language"))

5.4 Read and write

rw_bm <- brm(
  read_write ~
    language * dim_1,
  data = lang_use,
  family = bernoulli(),
  # Same priors as understand_bm
  prior = understand_priors,
  backend = "cmdstanr",
  cores = 4,
  file = "./data/rds/rw_bm"
)
conditional_effects(rw_bm, effects = "dim_1:language")

5.5 Mean competence

lang_use %>%
  ggplot(aes(mean_comp)) +
  geom_density() +
  geom_rug()

lang_use %>%
  ggplot(aes(language, mean_comp)) +
  geom_jitter(width = 0.2, alpha = 0.5, height = 0) +
  expand_limits(y = 0)

lang_use %>%
  ggplot(aes(language, mean_comp, colour = dim_1)) +
  geom_jitter(width = 0.2, alpha = 0.5, height = 0) +
  expand_limits(y = 0)

lang_use %>%
  ggplot(aes(dim_1, mean_comp)) +
  geom_point() +
  facet_grid(~ language)

mean_comp_bm <- brm(
  mean_comp ~
    language *
    dim_1,
  data = lang_use,
  family = zero_one_inflated_beta(),
  backend = "cmdstanr",
  cores = 4,
  file = "./data/rds/mean_comp_bm"
)

mean_comp_bm
##  Family: zero_one_inflated_beta 
##   Links: mu = logit; phi = identity; zoi = identity; coi = identity 
## Formula: mean_comp ~ language * dim_1 
##    Data: lang_use (Number of observations: 578) 
##   Draws: 4 chains, each with iter = 1000; warmup = 0; thin = 1;
##          total post-warmup draws = 4000
## 
## Population-Level Effects: 
##                         Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept                   0.34      0.04     0.25     0.43 1.00     5510
## languageesperanto           0.95      0.13     0.71     1.20 1.00     5105
## dim_1                       0.14      0.06     0.03     0.26 1.00     5129
## languageesperanto:dim_1    -0.13      0.16    -0.45     0.18 1.00     5092
##                         Tail_ESS
## Intercept                   2906
## languageesperanto           3080
## dim_1                       3163
## languageesperanto:dim_1     2950
## 
## Family Specific Parameters: 
##     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## phi     4.31      0.27     3.81     4.88 1.00     5900     2633
## zoi     0.22      0.02     0.19     0.25 1.00     6071     2944
## coi     0.99      0.01     0.97     1.00 1.00     4188     1861
## 
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
conditional_effects(mean_comp_bm, effects = "dim_1:language")

5.6 Dim-1 and age

lang_use %>%
  ggplot(aes(age, dim_1)) +
  geom_point() +
  geom_smooth(method = "lm")

5.7 Rural vs urban

5.7.1 Understand

understand_rur <- emil_rur_clean %>%
  count(understand, ru_ur) %>%
  group_by(ru_ur) %>%
  mutate(
    ru_ur_prop = n / sum(n)
  ) %>%
  ungroup() %>%
  mutate(
    understand = factor(understand, ordered = FALSE),
    ru_ur = factor(ru_ur, levels = c("rural", "urban"))
  )
get_prior(
  n ~
    understand * ru_ur,
  data = understand_rur,
  family = Beta()
)
##                 prior     class                       coef group resp dpar
##                (flat)         b                                           
##                (flat)         b                 ru_ururban                
##                (flat)         b            understand50D50                
##                (flat)         b understand50D50:ru_ururban                
##                (flat)         b               understandAL                
##                (flat)         b    understandAL:ru_ururban                
##                (flat)         b                understandG                
##                (flat)         b     understandG:ru_ururban                
##                (flat)         b               understandVG                
##                (flat)         b    understandVG:ru_ururban                
##  student_t(3, 0, 2.5) Intercept                                           
##     gamma(0.01, 0.01)       phi                                           
##  nlpar bound       source
##                   default
##              (vectorized)
##              (vectorized)
##              (vectorized)
##              (vectorized)
##              (vectorized)
##              (vectorized)
##              (vectorized)
##              (vectorized)
##              (vectorized)
##                   default
##                   default
understand_rur_priors <- c(
  prior(normal(0, 3), class = Intercept),
  prior(normal(0, 1), class = b)
)

understand_rur_bm <- brm(
  ru_ur_prop ~
    understand * ru_ur,
  data = understand_rur,
  family = Beta(),
  prior = understand_rur_priors,
  backend = "cmdstanr",
  cores = 4,
  file = "./data/rds/understand_rur_bm"
)

understand_rur_bm
##  Family: beta 
##   Links: mu = logit; phi = identity 
## Formula: ru_ur_prop ~ understand * ru_ur 
##    Data: understand_rur (Number of observations: 9) 
##   Draws: 4 chains, each with iter = 1000; warmup = 0; thin = 1;
##          total post-warmup draws = 4000
## 
## Population-Level Effects: 
##                            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept                     -1.82      0.60    -2.93    -0.60 1.00     2295
## understandAL                  -0.47      0.66    -1.79     0.82 1.00     3290
## understand50D50                0.06      0.66    -1.31     1.29 1.00     2970
## understandG                    1.07      0.71    -0.42     2.35 1.00     2182
## understandVG                   1.38      0.74    -0.17     2.68 1.00     2080
## ru_ururban                    -0.53      0.56    -1.55     0.62 1.00     2697
## understandAL:ru_ururban        0.32      0.81    -1.29     1.85 1.00     3007
## understand50D50:ru_ururban     0.10      0.78    -1.50     1.57 1.00     3564
## understandG:ru_ururban         0.48      0.74    -1.08     1.87 1.00     3124
## understandVG:ru_ururban        0.57      0.73    -0.88     1.96 1.00     3422
##                            Tail_ESS
## Intercept                      2570
## understandAL                   2503
## understand50D50                2824
## understandG                    2140
## understandVG                   2591
## ru_ururban                     2745
## understandAL:ru_ururban        2585
## understand50D50:ru_ururban     3041
## understandG:ru_ururban         2653
## understandVG:ru_ururban        2795
## 
## Family Specific Parameters: 
##     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## phi    14.59     12.89     2.32    48.69 1.00     1439     1408
## 
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
conditional_effects(understand_rur_bm, "understand:ru_ur")

5.7.2 Speak

speak_rur <- emil_rur_clean %>%
  count(speak, ru_ur) %>%
  group_by(ru_ur) %>%
  mutate(
    ru_ur_prop = n / sum(n)
  ) %>%
  ungroup() %>%
  mutate(
    speak = factor(speak, ordered = FALSE),
    ru_ur = factor(ru_ur, levels = c("rural", "urban"))
  )
speak_rur_bm <- brm(
  ru_ur_prop ~
    speak * ru_ur,
  data = speak_rur,
  family = Beta(),
  prior = understand_rur_priors,
  backend = "cmdstanr",
  cores = 4,
  file = "./data/rds/speak_rur_bm"
)

speak_rur_bm
##  Family: beta 
##   Links: mu = logit; phi = identity 
## Formula: ru_ur_prop ~ speak * ru_ur 
##    Data: speak_rur (Number of observations: 10) 
##   Draws: 4 chains, each with iter = 1000; warmup = 0; thin = 1;
##          total post-warmup draws = 4000
## 
## Population-Level Effects: 
##                       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept                -1.97      0.36    -2.59    -1.18 1.00     1541
## speakAL                   0.83      0.47    -0.26     1.61 1.00     1773
## speak50D50                0.11      0.45    -0.87     0.94 1.00     2409
## speakG                    1.05      0.47     0.01     1.82 1.00     1456
## speakVG                   0.60      0.48    -0.50     1.42 1.00     1639
## ru_ururban               -0.22      0.39    -1.00     0.56 1.00     2276
## speakAL:ru_ururban        0.54      0.53    -0.59     1.55 1.00     2737
## speak50D50:ru_ururban     0.52      0.58    -0.73     1.58 1.00     2131
## speakG:ru_ururban        -0.02      0.53    -1.08     1.07 1.00     2379
## speakVG:ru_ururban        0.02      0.56    -1.10     1.11 1.00     2700
##                       Tail_ESS
## Intercept                 1655
## speakAL                   1909
## speak50D50                2259
## speakG                    2106
## speakVG                   2279
## ru_ururban                2604
## speakAL:ru_ururban        2566
## speak50D50:ru_ururban     2111
## speakG:ru_ururban         2193
## speakVG:ru_ururban        2523
## 
## Family Specific Parameters: 
##     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## phi    71.40     70.74     7.49   250.05 1.00      619     1074
## 
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
conditional_effects(speak_rur_bm, "speak:ru_ur")

5.7.3 Read and write

rw_rur <- emil_rur_clean %>%
  count(read_write, ru_ur) %>%
  group_by(ru_ur) %>%
  mutate(
    ru_ur_prop = n / sum(n)
  ) %>%
  ungroup() %>%
  mutate(
    ru_ur = factor(ru_ur, levels = c("rural", "urban"))
  )
rw_rur_bm <- brm(
  ru_ur_prop ~
    read_write * ru_ur,
  data = rw_rur,
  family = Beta(),
  prior = understand_rur_priors,
  backend = "cmdstanr",
  cores = 4,
  file = "./data/rds/rw_rur_bm"
)

rw_rur_bm
##  Family: beta 
##   Links: mu = logit; phi = identity 
## Formula: ru_ur_prop ~ read_write * ru_ur 
##    Data: rw_rur (Number of observations: 4) 
##   Draws: 4 chains, each with iter = 1000; warmup = 0; thin = 1;
##          total post-warmup draws = 4000
## 
## Population-Level Effects: 
##                          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept                    0.49      0.49    -0.70     1.34 1.00     1987
## read_writeYes               -0.99      0.61    -1.96     0.47 1.00     1301
## ru_ururban                  -0.07      0.54    -1.11     1.08 1.00     2066
## read_writeYes:ru_ururban     0.14      0.65    -1.26     1.29 1.00     1616
##                          Tail_ESS
## Intercept                    1981
## read_writeYes                2029
## ru_ururban                   2112
## read_writeYes:ru_ururban     2179
## 
## Family Specific Parameters: 
##     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## phi    36.85     49.71     1.94   171.77 1.00      534      902
## 
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
conditional_effects(rw_rur_bm, "read_write:ru_ur")